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Abstract — This paper introduces Bayesian supervised and 
unsupervised segmentation algorithms aimed at oceanic segmen- 
tation of SAR images. The data term, i.e,, the density of the 
observed backscattered signal given the region, is modeled by 
a finite mixture of Gamma densities with a given predefined 
number of components. To estimate the parameters of the class 
conditional densities, a new expectation maximization algorithm 
was developed. The prior is a multi-level logistic Markov ran- 
dom field enforcing local continuity in a statistical sense. The 
smoothness parameter controlling the degree of homogeneity 
imposed on the scene is automatically estimated, by computing 
the evidence with loopy belief propagation; the classical coding 
and least squares fit methods are also considered. The maximum 
a posteriori segmentation is computed efficiently by means of 
recent graph-cut techniques, namely the a-Expansion algorithm 
that extends the methodology to an optional number of classes. 
The effectiveness of the proposed approaches is illustrated with 
simulated images and real ERS and Envisat scenes containing 
oil spills. 

Index Terms — Oceanic SAR images. Segmentation, Markov 
Random Fields, Energy minimization. Graph cuts. Mixture of 
Gammas, Oil spills. 



I. Introduction 

A wide number of oceanic phenomena become visible on 
SAR images as they have distinct scattering characteristics, 
namely the sea surface roughness and thus the normalized 
radar cross section (NRCS). Among these phenomena are 
gravity waves, convective cells, oceanic internal waves, current 
and coastal fronts, eddies, upwelling processes, ship wakes and 
oil pollution [1]. The automatic detection of these signatures is 
of outmost interest for a panoply of ocean monitoring systems, 
both for security, for commercial and for research applications. 
An example of an automated ocean feature detection scheme, 
able to detect fronts, ice edges and polar lows, is described in 
[2]. Another example is the Ocean Monitoring Workstation 
(OMW), developed by the company Satlantic. Usually, the 
first two steps of the processing chains of such systems are 
segmentation followed by classification. The Segmentation 
step computes a set of regions defining an image partition, 
where the features of each region, for example the gray levels, 
are similar in some sense. The classification focus then on each 
region attaching a label to it. For example, oil spill detection 
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based on SAR images has bee approached by many authors 
with the referred to classification-segmentation scheme(see for 
example [3]); the SAR image is first segmented and then the 
classification is focused on the regions with lower scattering; 
the classifier then computes a number of features from these 
regions including shape, moments, scale parameters, etc. [3] 
based on which a decision on whether the region corresponds 
to oil or look-alike is taken. We should refer, however. 

An application where segmentation of oceanic SAR im- 
ages is important is oil spill detection. In fact, many ap- 
proaches to this issue have been proposed in recent years. 
Most of them follow the described segmentation-classification 
structure, although other methods exist that do not, like for 
example kernel-based anomaly detectors [4]. A review of 
SAR segmentation techniques for oil spill detection can be 
found in [3]. The present work, which is an elaboration of 
our previous works [5], [6], [7], describes algorithms for 
the segmentation of dark signatures in oceanic SAR images, 
following a Bayesian approach. The adopted technique uses as 
data model a finite Gamma mixture, with a given predefined 
number of components, for modeling each class density. In 
fact, this density is well suited to filtered intensity SAR images 
as shown, e.g., in [8] and in [9]. By using a mixture, we aim at 
describing the continuous backscattering variability that may 
be observed in the SAR sea data. Moreover, a mixture is able 
to describe densities presenting more than one maximum, as 
it is the case of oceanic multi-look SAR images histograms. 
This cannot be achieved by none of the statistical models for 
SAR intensities proposed in the literature, like for example 
any of the Pearson distributions. 

In this work we propose two supervised and one unsu- 
pervised algorithm. The supervised algorithms demand an 
interaction with the user that manually selects a region con- 
taining pixels from the dark signature of interest and a region 
containing water pixels. These regions should be representative 
and can be made up of different not connected parts. The 
unsupervised algorithm is an improvement of the supervised 
ones and is completely automatic. 

To estimate the parameters of the class conditional densities, 
a new expectation maximization (EM) algorithm was devel- 
oped. Details are described in appendix. When segmenting 
small sub-scenes of the image, a simplified data model with 
only one Gamma function per class can be used. 

The prior used to impose local homogeneity is a multi- 
level logistic (MLL) model Markov random field (MRF) [10], 
with 2nd order neighborhood. To infer the prior smoothness 
parameter controling the degree of scene homogeneity, we 
develop an EM algorithm that uses loopy believe propagation 
(LBP) [11]. We have also exploited different classic estimation 



methods, namely the least squares fit (LSF) and the coding 
method (CD) (see [10] for details) for comparison. 

To infer the labels, we adopt the maximum a posterior 
(MAP) criterion, which we implement efficiently and exactly 
with graph-cut techniques [12]. 

Although the segmentation of oceanic dark patches is 
typically based on the assumption of two classes, we have 
generalized the problem to an optional number of classes. The 
underlying integer optimization problem is now attacked with 
the graph-cut based a-Expansion algorithm [13]. 

To evaluate the accuracy of the proposed algorithms, differ- 
ent simulations addressing both the referred Gamma model as 
well as intensity images corrupted with Gaussian noise have 
been carried out for error rate assessment. 

The algorithms have also been applied to real SAR images 
containing well documented oil spills. For doing so, the scenes 
are divided in tiles and segmented individually. 

A. Related Work 

Related approaches to the problem of oil spill segmentation 
are built on off-the-shelf segmentation algorithms such as the 
adaptive image thresholding and the hysteresis thresholding. 
Entropy methods based on the maximum descriptive length 
(MDL) and wavelet based approaches have also been pro- 
posed. Another recently proposed segmentation methodology 
applies Hidden Markov Chains (HMC) to a multiscale repre- 
sentation of the original image. Hereby the wavelet coefficients 
are statistical characterized by the Pearson system and by 
the the generalized Gaussian family [14]. When other SAR 
products are available, for example polarimetric data, other 
methods have been described in the literature, like constant 
false alarm rate filters [15]. 

An example of an elaborated adaptive thresholding tech- 
nique is provided by [16]. In this method, an image pyramid 
is created by averaging pixels in the original image. From 
the original image, the next level in the pyramid is created 
with half the pixel size of the original image. A threshold 
is then computed for each level based on local estimates of 
the roughness of the surrounding sea and on a look-up table 
containing experimental values obtained from a training data 
set. 

Hysteresis thresholding has been used as the base for 
detecting oil slicks in [17]. The method includes two steps: 
applying a so-called directional hysteresis thresholding (DHT) 
and performing the fusion of the DHT responses using a 
Bayesian operator. The MDL technique, which basically con- 
sists in applying information theory in order to find the 
image description which has the lowest complexity, has been 
applied in [18] to segment speckled SAR images, namely those 
containing oil slicks. This segmentation method describes the 
image as a polygonal grid and determines the number of 
regions and the location of the nodes that delimit the regions. 
The two-dimensional wavelet transform, used as a bandpass 
filter to separate processes at different scales, has also been 
adopted to oil slick detection in the framework of an algorithm 
for automated detection and tracking of mesoscale features 
from satellite imagery. [2]. 



B. Contributions 

We approach oil spill segmentation using a Bayesian frame- 
work and a multi-level Logistic (MLL) prior. Several methods 
in the same vein have been proposed since the seminal work of 
Geman and Geman [19], see e.g., [20]. Applications of these 
ideas in the segmentation of SAR images can be found, e.g., 
in [21], [22], [23]. 

The main contributions of this work to the state-of-art in oil 
spill segmentation are the following: 

• the development of an EM algorithm to estimate the 
parameters of a mixture of a pre-defined number of 
Gamma distributions, in order to model the intensities 
in a SAR image 

• the development of an EM algorithm using LBP to 
estimate the smoothness parameter in the MRF used as 
pior in our framework 

• the application of recent graph-cut techniques for solving 
the energy minimization problem that arises from the 
followed Bayesian methodology. 

• the design of a semisupervised algorithm for oil spill 
segmentation supported on the tools referred to above. 

C. Paper Organization 

The article is organized as follows: Section 1 introduces 
the problem, with references to related work and giving the 
main contributions of the present work; Section 2 overviews 
the Bayesian methodology that builds the base to the proposed 
algorithms. In addition, it briefly reviews the concept of the a- 
Expansion technique that has been implemented to generalize 
the methodology to an optional number of classes; Section 3 
describes the used parameter estimation techniques; Section 4 
describes the supervised segmentation algorithms; Section 5 
describes the unsupervided algorithm and Section 6 presents 
results of segmenting simulated and real images applying the 
algorithms proposed in Section 4 and 5; Finally Section 7 
contains the main conclusions and future work remarks. The 
article also includes an appendix where the referred to EM 
algorithm, developed to estimate the parameters of the class 
conditional densities of the Gamma mixture data model, is 
described. 

II. Problem Formulation 

A. Bayesian Approach 

let C := {l,...,c} be a set of c classes and V := 
{1,2, . . . , A/"} be the set of N pixels (sites) where mea- 
surements y := {^1,^2, . . . ,^Ar}, the SAR intensities, are 
available. A labeling x := {xi, X2, . . . xn} is a mapping from 
V to £, i.e., it assigns to each pixel p e V 3. label Xp G C. 
Any labeling x can be uniquely represented by a partition of 
image pixels P = {Pi\l G jC}, where Pi = {p e V\xp = 1} is 
the subset of pixels to which the label / have been assigned. 
Since there is an one-to-one correspondence between labelings 
X and partitions P, we use these notions interchangeably. By 
applying a segmentation algorithm to the image y, we get 
X := {xi, f 2 • • • Xn}, where Xi is the inferred label for pixel 
i e V. 



B. Observation Model 

In our problem formulation, we assume conditional inde- 
pendence of the measurements given the labels, i.e., 



N 



p{y\x) = Y[p{yi\xi) = n n ^(^^1^^)' 



(1) 
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where p(-|0O is the density corresponding to class / and (j)i 
the correspondent vector of parameters. The adopted density 
is a finite Gamma mixture given by 



K 



p{y^\^') = Y.<P(y^\^'s)^ 



(2) 



where K is the number of Gamma modes in the mixture, 
i indexes the pixel, and, for the class /, O^ is the vector of 
parameters of the Gamma mode s, and a[ is the a priori 
probability of mode s. We denote 6^ := {0[, . . . , 6>^), a^ := 
(a^i,...,Qf^), and^^ := {a\0^). 
Given that p{yi\0\) is Gamma distributed, we have then 



P{ViWs) 



r(ai) 



y°' 'exp(-Aiyi), t/i > 0, (3) 



where ^^ := (a^,A^). The mean and variance of a random 
variable with the density dSji is, respectively, a^/Ag and 

The mixture parameters 0^ are estimated from the data 
by applying the EM Gamma mixture estimation algorithm 
described in Appendix. The procedure is the same both for 
the supervised and for the unsupervised algorithms. 

We now make a brief comment on our choice of the Gamma 
mixture for modeling the observation densities. Under the 
assumption of fully developed speckle, the complex radar 
amplitude is zero-mean circular Gaussian distributed [24] . The 
average intensity computed over a number of independent 
random variables with the same density is, thus. Gamma 
distributed. It happens, however, that in sea SAR imaging 
one or more of the above assumptions may fail, rendering 
the Gamma density a poor model for SAR intensity [14]. A 
line of attack to obtain better models is to use more flexible 
parametric families, such as the Pearson System [25, Ch. 4.1] 
(see also [9]). However, in our problem each class density 
exhibits much more variability than that of accommodated 
by the the Pearson System. We have very often, for exam- 
ple, multi-modal densities. We should resort, therefore, to 
a mixture of Pearson System densities, what would lead to 
complex learning procedures. We have experimentally ob- 
served, however, that the Gamma mixture yields very good 
fittings for real SAR histograms, obtained with a moderately 
complex learning algorithm. For this reason, we have adopted 
the Gamma mixture model. 

C. Prior 

A second assumption we are making is of local smoothness 
of the labels in a statistical sense. It is more likely to have 
neighboring sites with the same label than the other way 
around. We model this local smoothness with a second order 



MRF, P {x)\ more specifically, by an MLL model (Ising model 
in the case of two classes). The Markov property assumes that 



p{xi\ XjJ e V) =p{xi\ XjJ eMi), 



(4) 



where Mi is the set of neighbors of site i. If 
p{xi\ Xj.,j ^ J\fi) > 0, then the Hammersley-Clifford 
Theorem states that p{x) has the form of the Gibbs 
distribution 

p(x) = iexp-^(-), (5) 

where Z is the so called partition function and U is the energy 
function 



U{x)= ^Vci{a 



(6) 
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where C is the set of cliques and Vd (x) is the clique potential 
defined over clique d. In this work, we have pair- wise cliques 
defined on a second-order neighborhood (8 pixels). That is, 
C = {{ij) : i e Mj, j eMi, i> j). 

In these conditions, the MLL clique potentials, in the 
isotropic case, is given by 



Vcl {Xr, Xs) = -PS {Xr - Xg) , 



(7) 



where r, 5 G cl, S{x) := Io{x) is the indicator function of 
set {0}, and parameter (3 > controls the degree of scene 
homogeneity. 



D. Maximum a Posteriori Estimate 

The posterior of the labeling given the observed data is 

p{x\y) (X p {y\x) p (x) . (8) 

In order to infer x, we adopt the MAP criterion. This amounts 
to maximize the posterior density of the labeling given the 
observed data: 



X = argmaxp {x\y) , 



(9) 



and is equivalent to minimize the negative logarithm of ([8]). In 
this sense, we may rewrite the problem in the following way: 



where 



X = argmin£^(xi, . . . ,Xn), (10) 



E{xu...,xn) = -logp(x|y)+c^^ (11) 



and c^^ denotes an irrelevant constant. From equation ([8]), we 
have 
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E(xi,...,x^) = ^E^(x,)+ ^ E'^^{x,,xj), (12) 



with 



E'{x,) = -\ogp{y,\x,) (13) 

E'^\xi,Xj) = -f36{xi-Xj). (14) 



E. Energy Minimization 

As already stated, we are concerned with the minimization 
of E {xi^. . .xn) given by ([T2b . which we term energy. 
For two classes, {i.e., c = 2), the global minimum of 
E {xi^. . . xn) can be computed exactly by applying the graph- 
cut algorithm described in [12]. This is a consequence of en- 
ergy being graph-representable, i.e., E^'^ (0, 0) + £^*'-^ (1, 1) < 
E^'^ (0, 1) + E^^^ (1,0) (for details see [12]). For more than 
two classes, the solution of (|9]) can be approximately computed 
by the a-Expansion technique [13], also based on graph- 
cut concepts. This algorithm finds the local minimum of the 
energy within a known factor of the global minimum. 

We now give a brief description of the a-Expansion al- 
gorithm. Given a label a, a move from a partition P (with 
correspondent labeling x) to a new partition P' (with corre- 
spondent labeling x') is called an a-Expansion if Poc C P'^ and 
P/ C Pi for any label / / a. In other words, an a-Expansion 
move allows any set of image pixels to change their label to 
a. The algorithm cycles through the labels in V in some fixed 
or random order and finds the lowest a-Expansion move from 
the current labeling. If the expansion move has lower energy 
than the current labeling, then it becomes the current labeling. 
The algorithm terminates when no a-Expansion move exist 
corresponding to a local minimum of the energy. 

III. Estimation of Parameter p 

In this work, we consider three different techniques to 
determine the smoothness parameter /3: a new EM method 
hereafter introduced that uses loopy belief propagation (LBP) 
and the classical LSF and CD methods. Because LSF and CD 
assume the existence of labeled data, we have conceived an 
iterative labeling-estimation scheme, which alternates between 
a labeling step and an estimation step until the convergence of 
P is attained. On the contrary, the EM estimation algorithm, 
that we have called "Loopy-/3-Estimation", is a one-shot 
technique. In the following sections, we briefly review the 
LSF and CD methods and provide a detailed description of 
the Loopy-/3-Estimation method. In this section, the class 
parameters (j) := (^^, ^^, . . . , ^^) are assumed known. 

A. Least Squares Fit 

This procedure for parameter estimation in MRF, described 
in detail in [10], is based in the following equation that holds 
for the MLL model, for every pixel p, with neighborhood ATp, 
and for every label pair Xp^x^ e C : 



/3 [n{xp) -n{xp)] 



log 



'p{xp\xxp,y) 



^p{x'p\xM^,y) ) 

- {E^ (x;) - EP {xp)) , (15) 

where n{xp) is the number of pixels in the set JVp with 
the same label as Xp and xj^f^ := {xi^ i G JVp}. We use 
histograms to estimate the joint probabilities p {xp\xj\fp^y) 
and p (^Xp\x^p^y) in ([T5]) : assuming that there are a total of 
M distinct 3x3 blocks in the image lattice with a given label 
configuration xj^f^, then we take 

y^ {xp\xj^p,y) 



where H {xp\xXp^y) is the number of times that a particular 
3x3 configuration {xp\xXpjy) occurs. The expression is then 
evaluated for a number of distinct combinations of Xp , x^ and 
xj\fp in order to obtain an over-determined linear system of 
equations that is solved in order to /3. 

B. Coding Method 

In this method the key idea (see [10]) is to partition the set 
V into sets V^^\ called codings, such that no two sites in one 
-pik) ^j,g neighbors. In the present work, the neighborhood ATp 
is of 2nd order thus yielding four codings. As the pixels in 
-pik) ^^Q. j^Q^ neighbors, the variables associated with these 
pixels, given the labels at all other pixels, are mutually 
independent under the Markovian assumption. The following 
simple product is thus obtained for the likelihood: 

p(^)(x|/3,y)= \[ p{xp\xMp,P,y), (17) 

pev^^^ 

Maximizing (fTTl) in order to P gives the coding estimate 
P^^\ Although it is not clear how to combine the results 
optimally, the arithmetic average, as suggested in [10], is an 
intuitive scheme that was adopted in this work. 

C. Loopy- j5 -Estimation 

We are seeking 13ml = argmax/3p(?/|/3), the ML estimate 
of the smoothness parameter /3. Based on the fact that the 
marginal density p{y\P), the so-called evidence, is a sum over 
the missing labels x, i.e., 



p{y\P) = ^p{y^x\P) 

X 

= ^p{y\x)p{x[ 



(18) 



we develop an EM algorithm [26] for the ML estimation 
of the parameter /3. The EM algorithm alternates between 
two steps: the E-step computes the conditional expectation 
of the logarithm of the complete a posteriori probability 
function, with respect to the missing variables, based on the 
actual parameter value; the M-step updates the value of the 
parameter, by maximizing the expression obtained in the 
E-step with respect to that parameter. We now derive the 
E-step and the M-step. 



p[xp\xM^,y) = 



M 



(16) 



E-step: 










Q{l3;l3t) = E[logp{y,x\l3)\y,l3t] (19) 


= E[logp{y\x)\y,/3t] (20) 


+E[\ogp{xmy,pt]. (21) 


Recalling that the MLL prior is given by 


pN/3)=2(^)exp 


p > ; s{xi-xj) 


, (22) 


with 


Z(/3) = )_Jexp 


p >; sixi-x,) 


(23) 
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^ij(Xi,Xj) 
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Fig. [T] depicts the graph that represents our pairwise MRF 
for computing the two-node behefs. In the square lattice, 
ipij {xi , Xj ) stands for the interaction potential that penalizes 
every dissimilar pair of neighboring labels, and (j){xi^yi) = 
p{yi\ Xi) represents the statistical dependency between the 
labels Xi and the measurements yi. For computing hij (xi^Xj), 
(j){xi,yi), is set to a constant value, making the result inde- 
pendent of the data values y. 

We solve ([28l) by a line search type algorithm, ensuring that 



< 0, thus corresponding to a maximum of Q. 



Fig. 1. Lattice representing the pairwise MRF. 



we obtain, up to an irrelevant constant. 



Q(/3;A) = -log^(/3) + /3 Y. E\b{x,-x,)\y,^t\. 



M-step: 



(3t+i = arg max Q (/?, A) . 



(24) 
(25) 



The stationary points of Q are the solution of 

dQ _ d\og[Z{/3)] 
dp dp 



(26) 



+ X! X! Hxi-Xj)p{xi,Xj\ y,/3t))=0. 



By introducing expression (1231) into (1261) and if we consider 
that 5 {xi — Xj) takes non-zero values only for equal labels, 
we obtain 



dQ 

dp 



i,jGcl k—1 

-p{xi = k,Xj = k\ Pt) = 0. 



(27) 



Since computing exact marginal distributions is infeasible 
in our case, we replace them by pseudo-marginals using BR 
This approach has been successfully applied in problems of 
approximate parameter learning in discriminative fields [27]. 
BP is an efficient iterative algorithm in which local messages 
are passed in graphical models. For singly-connected (loop 
free) pairwise MRFs, the two-node beliefs will correspond 
to the exact two-node marginal probabilities. In our case, 
however, the graph that corresponds to the MRF contains 
cycles, preventing the basic BP algorithm to be applied. 
We resort to a slightly modified version called loopy belief 
propagation. In practice, this algorithm has often delivered 
good results [11]. We approximate the marginal probabilities 
with the two-node beliefs bij (xi^Xj) and b^j {xi^Xj). These 
will provide approximations respectively for the marginals 
p{xi^Xj\l3t) and p{xi^Xj\y^ Pt). By doing so, the M-Step of 
our EM algorithm is given by the sum of the differences 
between the two node beliefs that take the evidence y into 
consideration and those that do not make use of y\ 



^ i,jed fe=l 



■ bij {k, k) = 0. 



(28) 



D. A few Remarks about the Vector of Parametres (j) 

In the previous section, we have considered that the class pa- 
rameter vector (j) is known and only the smoothness parameter 
P is to be inferred. However, still using the beliefs computed 
by LBP, the vector (j) could have been inferred simultaneously 
with p by including in the Q function the additive term (l2Qb 
corresponding to the class densities. We would have obtained 

Q{^',^uPt) = E[\ogp{y\x)\y,(l)upt] 

N c 

= X]X]logp(?/i|^^)p(xi =l\y,(t)t,pt), 

where the probabilities p{xi = l\y^(j)t^Pt) are given by the 
LPB method. The M-step would consist then in two decoupled 
maximizations; one with respect to /3 and another with respect 
to (j). This approach is, however, beyond the scope of this 
paper. Nevertheless, we present in Section |Vl an unsupervised 
algorithm, which is suboptimal but faster than the LBP based 
EM algorithm aimed at the inference of both the parameters 
P and (j). 

An alternative to the proposed EM scheme based on LPB 
is using Monte Carlo techniques to cope with the difficulty 
in computing the E-step and the partition function Z [28], 
[29], [30]. Supported on the performance of the LBP, we 
believe, however, that the EM scheme based on LPB is, for 
the present problem, much more faster than the Monte Carlo 
based techniques. 

IV. Supervised Segmentation 

In this section, we introduce two supervised algorithms 
aimed at the segmentation of sea SAR images. In both algo- 
rithms, the first step is the estimation of the class parameters 
used in the data model. This is done by asking the user to de- 
fine representative regions of interest (ROI) in the image. Once 
the regions are defined, different approaches may be followed: 
if a Gamma mixture is assumed for the observed data, as in 
the case of SAR images, the EM class parameters estimation 
algorithm described in Appendix is applied to infer the class 
conditional densities. If we are segmenting small sub- scenes 
and the radar range spreading loss has been compensated in the 
underlying SAR images, one single Gamma function per class 
often provides a good modeling for the SAR intensities. In this 
case, a common ML Gamma estimator is used instead of the 
EM procedure. After this step, the data model is considered 
to be known and is used thereafter. 



The pseudo-code for two supervised algorithms is presented 
in Algorithm 1 and Algorithm 2. Algorithm 1 is of generalized 
likelihood type [31] implementing an iterative labeling scheme 
with two steps being performed alternately: the /3 -Estimation 
and the segmentation. Algorithm 2 is a one-shot procedure that 
performs /3-Estimation using the Loopy /3-Estimation method 
described in Section IIII-CI 

Algorithm 1 Supervised Segmentation Using LSF/CD (3- 

Estimation 

Require: Initial parameter (3 = Pq and estimated class pa- 
rameters (j) 
Compute eLabel = E^{xp) for every pixel p using ^ 



while 



A/3 



< 5 OX nrlterations < NrlterationsMax do 



Compute X 
Compute P 

end while 

return (x,/3 



a-Expansion(;^, eLabel) 
l3-Estim3.tion{x^ eLabel) {LSF or CD}. 



Algorithm 2 Supervised Segmentation Using Loopy-/3- 

Estimation 

Require: estimated class parameters cj) 
1: Compute eLabel = E^{xp) for every pixel p using (p 
2: Compute /3 = Loopy-/3-Estimation(eLa6eO 
3: Compute X = GraphCut-Segmentation(/3, eLabel) 
4: return (x,^) 



V. Unsupervised Segmentation 

The unsupervised method is an improvement of the super- 
vised algorithms described in the previous sections. The main 
difference is that the data model is not considered to be known 
but is also iteratively estimated along with the smoothness 
parameter and the segmentation. The scheme needs a rough 
initialization of the data model parameters. We computed this 
initialization by fitting an EM Gamma mixture to the complete 
data. 

As our EM algorithm automatically eliminates unnecessary 
modes, we start with an overestimate of K, the number of 
modes in the mixtures. In our experiments with real data, the 
maximum number of modes we got was four. 

Then, different strategies are possible. In the experiments 
reported in the next section, when considering c = 2 classes, 
we have assigned the mode with a lower mean value to one 
of the classes and the remaining modes to the other class. 

In each iteration, we compute the ML estimate of the 
vector (j) based on the previous segmentation, compute the new 
segmentation, and finally the new value of f3. The algorithm 
may use any of the three /3— parameter estimation methods and 
is applicable to an optional number of classes. The algorithm 
stops when both the parameter (3 and the class parameters 
converge to stable values. 



Algorithm 3 Unsupervised Segmentation 

Require: arbitrary parameter /3 = /3o, initial class parameters 
(j) = (j)Q (EM Gamma Mixture Estimation) 
1: Compute eLabelo = EP{xp) for every pixel p, using 0o 
2: Compute initial labeling 

X = xo = a-Expansion(/3o, eLabelo) 
3: for Stop criterium is not met do 
4: Compute (/) =ML-Estimation(f) 
5: Compute eLabel = E^{xp) for every pixel p, using ^ 
6: Compute X = Qf-Expansion(/3, eLabel) 
7: Compute $ = /3-Estimation(x, eLabel) {use LSF, CD 

or Loopy}. 
8: end for 
9: return (x,/3,^j 



VI. Results 

This section presents results of applying the proposed 
methodology to simulated and real SAR images, as well 
as some considerations regarding time complexity of the 
algorithms. 

A. Simulations 

We have performed simulations corresponding to Gamma 
data terms with two classes, and evaluated the overall accura- 
cies of the obtained segmentations using Algorithms 1 , 2 and 
3. We also illustrate in detail the EM Gamma mixture esti- 
mation algorithm, by applying it to an example of simulated 
data. 

1) Segmentation Results with Gamma Data Term: In order 
to compare the performance of the three proposed methods, 
we have tested Algorithms 1, 2 and 3 on simulated images 
generated by adding Gamma noise to ground-truth images 
containing two classes. We have used the ground-truth de- 
picted in Fig.|2](a) and added noise with Gamma distributions 
having mean values of five and nine. The parameters of the 
distributions were choosen in order to obtain increasing values 
of variance a'^, corresponding to noisier images. We have then 
applied Algorithm l,with LSF and CD methods. Algorithm 2, 
and Algorithm 3 with Loopy-/3-Estimation. Furthermore, for 
comparison, the images were also segmented tuning the beta 
value manually and for the case that no prior is used (13 = 0). 
Fig. [2] shows the segmentations obtained for a = 2.6. As we 
can see from /3 = this is a hard problem, on which LSF and 
CD fails the inner structures, but Loopy-/3-Estimation provides 
better results, both in the supervised as in the unsupervised 
method. This behavior is further confirmed by the rank in Fig. 
[51 which gives the OA obtained for images with different a 
values (corresponding to images with more noise) for the six 
different segmentation processes referred above. 

We have also tested the proposed algorithms in images 
generated by adding Gamma noise to a ground-truth similar 
to an oil patch (see Fig. \^. Two images. Image A and Image 

B, with the histograms depicted in Fig. IH corresponding to 
different segmentation difficulty levels are segmented. The best 
segmentation possible in this framework, achieved by tuning 



GroundTruth 



Segmentation for best 
beta tuned manually 




(a) 

Segmentation 
using 'Algorithm-2' 



(b) 

Segmentation for 

beta = 




Segmentation using 
'Algortihm-r with LSF 



(d) 

Segmentation using 
Algortihm-r with CD 




Segmentation using 
'Algortihm-3' with Loopy Estimation 




Fig. 2. Ground-truth and results of different segmentation processes for an 
image with Gamma noise having mean values of five and nine and with a a 
value = 2.6. Notice the good performance of Algorithm 2 and Algorithm 3, 
implementing the Loopy-/3-Estimation. 



the /3 value manually, and the segmentation obtained with 
no prior information (setting /3 = 0) are also displayed for 
comparison. Fig. [5] shows the results obtained for Image A 
and Fig. [6] the results obtained for Image B. For Image A, the 
segmented images using Algorithm 1 , with CD and with LSF, 
are not shown, as they are almost equal to the image segmented 
with Algorithm 2. For Image B, Algorithm 3 has not provided 
good results and the segmentation is not displayed. The bad 
performance of Algorithm 3 in this case arises from the not so 
good initialisation of the class parameters, due to the complete 
overlapping from the oil and the water histogramms. 

2) EM Algorithm for Gamma Mixture: In this subsection 
we illustrate the behavior of the EM algorithm designed to 
infer Gamma mixtures. For details on the theoretical issues, 
we refer to the Appendix. In Fig. [71 we can see the ground- 




Fig. 3. OA for images corrupted with Gamma noise with increasing a values: 
TM = for /3 tuned manually; LE = Algorithm 2, using Loopy Estimation; 
LSF = Algorithm 1 using LSF; CD = Algorithm 1 using CD, NP = no prior 
and UNS = Algorithm 3 using Loopy Estimation. 
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Fig. 4. Gamma Model: (a) Ground-truth, (b) Histogram of Image A, 
(c) Histogram of Image B, (d) Image A, (e) Image B with superimposed 
delimiting line, for better visualization 
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Fig. 5. Segmentation results for image A: (a) No prior information, (b) f3 
value tuned manually, (c) Segmentation using Algorithm 3 with Loopy-^- 
Estimation, (d) Segmentation using Algorithm 2 (Algorithm 1 provided the 
same results with only 0.1% difference in OA), (e) graph depicting the ^ 
estimation in Algorithm 2. 



truth used for simulating the image and the generated image, 
according to the densities shown in Fig.[8l The two classes, oil 
and water have been modeled respectively by a mixture of two 
Gamma functions and a mixture of three Gamma functions. 
For this particular case, we have selected a ROI containing 
178 pixels for representing the water and a ROI containing 142 
pixels for representing the oil. After 20 iterations of the EM 
algorithm, we obtained the approximations for the probability 
distributions depicted in Fig. [9] and Fig. \TU\ The obtained 
results are quite reasonable for the small sample sizes used. 

B. Real Images 

We have tested the proposed methodology with three real 
SAR images for the special purpose of oil spill detection. For 
this type of application, an unsupervised algorithm is more 
indicated, and so we decided to apply Algorithm 3. Because 
Loopy Estimation has proven to be an effective method in 
the simulations, we have choosen Algorithm 3 with Loopy 
Estimation. Nevertheless, in order to have a comparision 
unsupervised approach versus supervised approach, we also 
applied Algorithm 2 to two of the images. In oil spill detection. 
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Fig. 6. Segmentation results for Image B: (a) No prior information, (b) (3 
value tuned manually, (c) graph depicting (3 estimation using Algorithm 2, (d) 
Segmentation using Algorithm 2, (e) Segmentation using Algorithm 1 with 
LSF and (f) Segmentation using Algorithm 1 with CD. 





(a) 



Fig. 7. (a) Ground-truth used for simulating an oil spill. Black represents 
oil and white water; (b)Simulated SAR image 



the number of classes is typically set to c = 2, although more 
classes can be considered if we are interested in distinguishing 
other phenomena occuring at the same time in the region of 
interest. The application of the algorithms is straightforward: 
the key idea, like in most state-of-the-art oil spill detection 
methods (see for example [16]), is to partition the image in 
tiles and run the algorithm separately for each part. This step 
is preceeded by the application of a landmask to the image, 
what can be done using external coastline information or by 
adopting some coastline self-extraction procedure. After the 
segmentation of each tile, a procedure for grouping patches 
detected on the tile borders should be carried on. Another 
possibility to increase segmentation coherency in the borders 
is to define overlapping tiles or to force continuity to some 
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Fig. 8. Probability functions used to generate the simulated image with 
superimposed histogram of generated data set. A three modes function for 
water and two modes function for oil was used. 
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Fig. 9. True and estimated class densities for oil. 
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Fig. 10. True and estimated class densities for water. 




Fig. 11. ERS-1 image from the Sicily Channel, Italy, acquired on 30th 
January of 1992. The smaller and larger squares are sub-scenes that have been 
used respectively to estimate the wind direction and to apply our algorithm. 
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Fig. 12. Closer look to the sub-scenes of the ERS-1 image. Part (a) is used 
for wind direction estimation; the estimated direction has been overlapped as a 
green arrow. Part (b) contains an oil spill. The coloured regions correspond to 
the ROI's selected for class parameter estimation in the supervised algorithm 



degree on the estimated class and/or smoothness parameters 
from one tile to the next. Nevertheless, at the moment we 
are not doing this and these are considered future possible 
improvements to our methodology. We segmented oil spills 
contained in three different scenes, described in the following 
subsections. 

1) Segmentation of an ERS-1 sub-scene: We have seg- 
mented part of an ERS-1 image from the Sicily Channel, 
Italy, that has been acquired on the 30* January of 1992. 
This image is referred in the ESA web pages regarding oil 
slicks http : //earth . esa . int /ew/oil_slicks/ and 
contains three oil slicks, along with information regarding 
wind direction and intensity and existence of ships, ship wakes, 
natural oil films and currents. 

Fig. [TT] provides a quicklook of the scene with two 
squares superimposed: the larger representing the part to be 
segmented and the smaller representing a part used for wind 
estiamtion. Fig. [121 provides zooms of the referred squares. 
By applying the Radon Transform to the smaller square, the 
estimated direction has been calculated and is depicted in the 
image. The direction was consistent with the measured value 
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Fig. 13. Fitting of a mixture of Gammas to the data from the red square in 
the ERS-1 image. 




Fig. 14. Segmentation of ERS-1 subscene containing oil with (3 
estimated using Algorithm-3 with the Loopy-^ -Estimation method. 
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reported in the url site and is at the origin of the well-known 
"feathering" effect that can be observed in this linear spill. 
We have computed the backscattering values of the image, by 
performing calibration using the ESA provided BEST software 
Chttp : //earth. esa. int /services /best/) and 
then applied Algorithm-3 with Loopy-/3-Estimation. The 
result of the Gamma mixture estimation, in the initialisation 
of the algorithm, is depicted in Fig. [131 After only three 
iterations, both the class parameters and the smoothness 
parameter have converged. A /3 value equal to 1.4 was 
estimated. The segmentation is displayed in Fig. [141 We also 
applied Algorithm-2: we selected two ROI's (160 pixels for 
water and 77 for oil) in the image (shown in Fig. [12]) and 
computed ML Gamma Estimators for the two classes (see 
Fig. [15]). We then applied the algorithm and obtained an 
estimated P value equal to 1.44. The segmentation result is 
given in Fig. [16] and is similar to the result obtained with 
Algorithm-3. 

For comparison, we also provide the segmentation obtained 
with no prior (corresponding to /3 = 0) in Fig. [TT] In practice, 
the estimated value has proved to deliver a good segmentation. 
When lower /3 values were used, the result was noisier and 
for higher /3 values, the details of the spill disappeared. 
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Fig. 15. Class Parameters Estimation for Algorithm-2. 




Fig. 16. Segmentation of ERS-1 subscene containing oil with ^ 
estimated using Algorithm-2. 
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Fig. 17. Segmentation of ERS-1 subscene containing oil, with /3 = 0. 




Fig. 18. (a) ASAR image fragment with coloured regions for class parameter 
estimation in Algorithm-2; (b) Segmentation result applying Algorithm-2; (c) 
Segmentation result applying Algorithm-3 with Loopy-/3-Estimation. 




Fig. 19. Display of the ASA WSM full resolution image of the Prestige 
accident occured in November 2002 in Galicia. 



2) Segmentation of an Envisat ASAR IM sub-scene: 
We have also applied Algorithm 3 to a fragment of an 
ASAR Image Mode image, acquired on 19 July 2004, in 
the ocean between Cyprus and Lebanon. The fragment 
contained an occurred oil spill of circa 10 km's (centered 
on ^ 33°N,33°E39') that was documented on the EC 
Oceanides project (a project in the framework of the Europe 
's Global Monitoring for Environment and Security initiative) 
database. After six iterations, we achieved convergence of 
the P and values. The segmentation result is given in Fig. 
[TSl corresponding to a estimated /3 = 1.83. The unsupervised 
segmentation result is compared with the one provided by 
Algorithm-2, where the user provided ROI's for water and 
for oil that were used to estimate the class parameters at 
the beginning of the process (see Fig. [18]). In this case a 
/3 = 1.75 was obtained. 

3) Segmentation of an Envisat ASAR WSM sub-scene: In 
order to demonstrate the viability of applying Algorithm 3 to a 
whole ASAR_WSM scene, we have run it over the very well 
known image of a confirmed oil spill, namely the Prestige 
case (see Fig. [19]). This accident took place in November 
2002, in Galicia (Spain), when a tanker carrying more than 
20 million gallons (around 67,000 tons) of oil split in half 
off the northwest coast of Spain on 19 November 2002, 
threatening one of the worst environmental disasters in history. 
For segmenting the image, we have first partitioned it in tiles of 
600x600 pixels each and then applied Algorithm-3 to each tile 
independently. No post-processing of the borders or "clean-up" 
operations (like for example morphological operations) have 
been carried out. Fig. [20] shows the segmented image, made 
up of the concatenation of the segmentation of the individual 
tiles. As we can see, the result can be considered in general 
very good, with only some discrepancies located on the tiles' 
borders. For initialising the algorithm, the EM Gamma mixture 



estimation procedure was applied with four modes {K = 4). 
In fact this number is enough to provide a good fitting of the 
intensity data of the SAR image. On the other hand, our EM 
algorithm allows us to start with a higher mode number, and 
decreases this number automatically. We depict the results of 
this fitting on one of the tiles, shown in Figure [23] by showing 
the histogram of the tile with the superimposed estimated 
Gamma mixture (see Fig. [21] (a)). As we have explained, the 
mode corresponding to the lower mean value is assigned to 
the oil class and the others to the water class, providing an 
initialisation to the class parameters. The estimations of these 
parameters are actualised along the algorithm and, after nine 
iterations, we obtain the distributions shown in Fig. [21] (b). 

To fully demonstrate the possibilities of our algorithm, 
we again run it on the tile shown in Figure [23] but this 
time setting the number of classes to c = 3. By doing so, 
we hope to be able to segment a third ambiguous zone, 
corresponding to intermediate radiometry levels and probably 
due to atmospheric conditions originating a front. In this case 
we choose to fitt a Gamma mixture of c modes to the data 
for initialising the class parameters, as depicted in Fig. [22] 
With this initialisation, the first obtained segmentation using 
Algorithm 3 is displayed in Fig. [24] after only nine iterations, 
all parameters have already converged to the segmentation 
given in Fig. [25] When comparing this segmentation with 
the one provided by a state-of-the-art algorithm, namely by 
multiscale HMC model in [14], we consider to have obtained 
a very good result. 

4) Time considerations: The proposed algorithms display, 
in practical cases, a computational complexity of 0{N), with 
N being the number of pixels in the image. In fact, the most 
time-consuming steps in the algorithms are the a-Expansion 
and Graph-Cut segmentation routines. These use the max- 
flow code implementation referred in [32] that has complexity 
0{n'^-^) in the worst-case, however, in most situations is 
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Fig. 20. Display of the segmentation results of applying Algorithm 3 with 
Loopy-/3-Estimation to the image displayed in Fig. [19] 





Fig. 22. Histogramm of the data values (normalised) with superimposed 
estimated Gamma mixture when the number of classes is set to c = 3 




Fig. 23. Display of one of the individual tiles in which the image in Figure 
[191 has been divided for segmentation. 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



(b) 



Fig. 21. Results of applying the EM Gamma mixture estimation procedure to 
the tile depicted in Figure [23] (a) Histogramm of the data values (normalised) 
with superimposed estimated Gamma mixture; (b) Distributions corresponding 
to the oil class (left) and water class (right) after nine iterations of Algorithm 
3: Histogramm and superimposed fitting. 




Fig. 24. Display of the initial segmentation result of applying Algorithm 3 
with Loopy-^-Estimation and three classes to the tile depicted in Figure [23] 




Fig. 25. Display of the final segmentation result, after nine iterations, 
of applying Algorithm 3 with Loopy-/3-Estimation and three classes to the 
bottom left tile of Fig. 17. 



0{N). Both Algorithm 1 and Algorithm 3 converge in few 
iterations and, as an example, for images of the dimension used 
in the first simulation (64x64) (see lVI-ATTl) , only a few seconds 
are needed to run the algorithms. For a 600x600 pixel tile, like 
those in the Prestige example in lVI-B.3l Algorithm 3 is taking 
^ 1 minute, but this time performance could be improved by 
a more efficient implementation of our procedure, namely by 
improving the implementation of the Loopy-/3-Estimation and 
by migrating our matlab code fully to C-code. 

VII. Conclusions 

The results of applying the proposed methodology to simu- 
lated images with Gamma data models and to real SAR data 
are promising. The developed EM Gamma mixture estimation 
algorithm, when incorporated into the proposed segmentation 
algorithms, has proved to be an efficient tool for data modeling 
in SAR intensity images. 

With the supervised algorithms, high OA accuracies have 
been achieved even for simulated images with high levels 
of noise. In general the Bayesian approach has resulted in 
an OA increase in the segmentation process between 10 and 
15%, when compared to the segmentation using no prior 
information. Hereby, Algorithm 1 with LSF and Algorithm 1 
with CD methods provided similar results, with a performance 
close to that obtained by setting the /3 value manually. On the 
other hand, the LSF estimation procedure seems to become 
less reliable for noisier images. When compared to Algorithm 
2, using Loopy-/3-Estimation, we see that this last method 
allways provides equal or better results, outcoming Algorithm 
1 . Another advantage of Algorithm 2 is that it is usually faster 
than Algorithm 1, being a one-shot process. Both algorithms, 
by introducing prior information into the segmentation pro- 
cess, increase the OA significantly. As a conclusion. Algorithm 
2 should be preferred to Algorithm 1 , when using a supervised 
method for oceanic sar images segmentation. 

Algorithm 3, totally unsupervised, has been conceived as 
an improvement to Algorithm 1 and 2. When testing it on 
simulated data the obtained results were very good, although 



slightly worst than the supervised ones, as expected. In gen- 
eral, after a few iterations, the class and smoothness param- 
eters converge to stable and meaningful values. By applying 
Algorithm 3 to real images containing documented oil spills, 
the segmentation has been considered successful! We could 
segment both linear and patch type oil spills. Furthermore, the 
applicability of the method to segment whole scenes, as well as 
to segment more than two classes, has also been demonstrated 
in a well known image from the Prestige accident. As a 
conclusion, we believe the presented methods are suitable to 
be used for segmenting oceanic SAR images. In particular for 
oil spill detection. Algorithm 3 seems to be a suitable method. 
Based on observations that oil spills in the ocean are often 
dragged by the wind and align more or less perpendicular to 
its direction, an interesting future issue is the incorporation 
of wind information into our segmentation algorithms. By 
adopting anisotropic MRF in the prior, we intend to reflect this 
directional dependency of the clique potentials. In the practice, 
the smoothness parameter (3 is no longer considered to be a 
constant value but is clique-dependent according to the wind 
direction and/or velocity. For estimating the wind profile from 
the SAR data, state-of-the-art algorithms are used. We have run 
promissing simulations, are currently testing this extension on 
real data, and expect to obtain interesting results on a near 
future. 

Appendix 
Fitting a Mixture of Gamma Densities 

As stated in the main text of this work, when segmenting 
real SAR images, the adopted data model for each class is 
a finite Gamma mixture. For completely defining the data 
model for each class (for lightness of the notation we have 
dropped the class index), we need to estimate 2K Gamma 
parameters, 6 = (^i, . . . ,6>x), with 6s = (as,As), and K 
a priori probabilities a = (ai, . . . , a^)). We infer (j) = 
(a^O) by computing its ML estimate from a training set. 
The ML estimate is computed via an EM algorithm [26]. 
Fitting a Gamma mixture is addressed in [33]. However, 
the authors consider TV-look SAR images meaning that the 
underlying random variables are the average of N independent 
and identically distributed exponential random variables, thus, 
having a Gamma density but with just one parameter free; if 
the mean is /i, then the variance is given by ji^ /N. We estimate 
both the mean and the variance for each Gamma distribution 
in the mixture, rendering the algorithm more adaptable to real 
measurements. 

The key point in the EM technique is the introduction of the 
so called missing data z, such that p{y\(t)) = J p {y ^ z\(j)) dz 
and p{y^z\(j)) is easier to manipulate than p{y\(j)). In the 
particular case of a mixture of densities, we will use as 
missing data a random variable, Zi, per site, with distribution 
p{zi = s) = as. It is interpretable as the probability of 
the s — th Gamma mode is selected at pixel i. The EM 
algorithm alternates between two steps: the E-step computes 
the conditional expectation of the logarithm of the complete 
a posteriori probability function, with respect to the missing 
variables, based on the actual parameter value. The M-step 
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updates the values of the parameters, by maximizing the ex- 
pression obtained in the E-step with respect to each parameter 
on turn, i.e., 



E-step 
M-step 



Q{^;^')=E{\ogp{y,z\^)\y,^'} (29) 
6^+^ =argmaxQ((^;(^^) (30) 



Denoting 



P{zi = s\yi,^* 



(31) 



and taking into account that XI <^s - 1' then (/)^+^ can be found 
among the stationary points of the Lagrangean 
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1 , (32) 



= 1 s=l 



+ log(7/,"^-^) - A,7/, + log(a,), (33) 

where A denotes a Lagrange multiplier. The expression for wl^ 
(see [34]) is given by 
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(34) 



In the M-step, after differentiating C in order to the unknown 
parameters and setting the derivatives to zero, we obtain a 
closed solution for the updating of the a priori probabilities 
ai's, but numerical iteration is needed for determining param- 
eters a^'s and A^'s of the Gamma densities. Expression ([35]) 
gives the update expression for a^'s. 
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(35) 



Equation (l36l) and (|37]) give the update expressions for the 
parameters A^'s and a^'s. 



A*+i = 






4+1 = ^-1 

where 



log iK) Ei=i wji + Ei=i log JyiWs 

Tias)' 



*K) 



(36) 

1 

(37) 
(38) 



is the psi function. We refer to the Appendix B of [35] for a 
very fast Newton procedure to compute the inverse of the psi(-) 
function. Expressions (l36l) and (|37]) are iteratively recomputed 
until convergence is obtained, starting from initial values 
computed from the observed data y. The initial parameter 
values are calculated in such a way, that the initial probability 
function is a sum of equidistant Gammas that span the most 
representative data range. The EM scheme converges in a few 
tens of iterations. 
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Abstract — This paper introduces Bayesian supervised and 
unsupervised segmentation algorithms aimed at oceanic segmen- 
tation of SAR images. The data term, i.e,, the density of the 
observed backscattered signal given the region, is modeled by 
a finite mixture of Gamma densities with a given predefined 
number of components. To estimate the parameters of the class 
conditional densities, a new expectation maximization algorithm 
was developed. The prior is a multi-level logistic Markov ran- 
dom field enforcing local continuity in a statistical sense. The 
smoothness parameter controlling the degree of homogeneity 
imposed on the scene is automatically estimated, by computing 
the evidence with loopy belief propagation; the classical coding 
and least squares fit methods are also considered. The maximum 
a posteriori segmentation is computed efficiently by means of 
recent graph-cut techniques, namely the a-Expansion algorithm 
that extends the methodology to an optional number of classes. 
The effectiveness of the proposed approaches is illustrated with 
simulated images and real ERS and Envisat scenes containing 
oil spills. 

Index Terms — Oceanic SAR images. Segmentation, Markov 
Random Fields, Energy minimization. Graph cuts. Mixture of 
Gammas, Oil spills. 



I. Introduction 

A wide number of oceanic phenomena become visible on 
SAR images as they have distinct scattering characteristics, 
namely the sea surface roughness and thus the normalized 
radar cross section (NRCS). Among these phenomena are 
gravity waves, convective cells, oceanic internal waves, current 
and coastal fronts, eddies, upwelling processes, ship wakes and 
oil pollution [?]. The automatic detection of these signatures is 
of outmost interest for a panoply of ocean monitoring systems, 
both for security, for commercial and for research applications. 
An example of an automated ocean feature detection scheme, 
able to detect fronts, ice edges and polar lows, is described in 
[?]. Another example is the Ocean Monitoring Workstation 
(OMW), developed by the company Satlantic. Usually, the 
first two steps of the processing chains of such systems are 
segmentation followed by classification. The Segmentation 
step computes a set of regions defining an image partition, 
where the features of each region, for example the gray levels, 
are similar in some sense. The classification focus then on each 
region attaching a label to it. For example, oil spill detection 
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based on SAR images has bee approached by many authors 
with the referred to classification-segmentation scheme(see for 
example [?]); the SAR image is first segmented and then the 
classification is focused on the regions with lower scattering; 
the classifier then computes a number of features from these 
regions including shape, moments, scale parameters, etc. [?] 
based on which a decision on whether the region corresponds 
to oil or look-alike is taken. We should refer, however. 

An application where segmentation of oceanic SAR im- 
ages is important is oil spill detection. In fact, many ap- 
proaches to this issue have been proposed in recent years. 
Most of them follow the described segmentation-classification 
structure, although other methods exist that do not, like for 
example kernel-based anomaly detectors [?]. A review of 
SAR segmentation techniques for oil spill detection can be 
found in [?]. The present work, which is an elaboration of 
our previous works [?], [?], [?], describes algorithms for 
the segmentation of dark signatures in oceanic SAR images, 
following a Bayesian approach. The adopted technique uses as 
data model a finite Gamma mixture, with a given predefined 
number of components, for modeling each class density. In 
fact, this density is well suited to filtered intensity SAR images 
as shown, e.g., in [?] and in [?]. By using a mixture, we aim at 
describing the continuous backscattering variability that may 
be observed in the SAR sea data. Moreover, a mixture is able 
to describe densities presenting more than one maximum, as 
it is the case of oceanic multi-look SAR images histograms. 
This cannot be achieved by none of the statistical models for 
SAR intensities proposed in the literature, like for example 
any of the Pearson distributions. 

In this work we propose two supervised and one unsu- 
pervised algorithm. The supervised algorithms demand an 
interaction with the user that manually selects a region con- 
taining pixels from the dark signature of interest and a region 
containing water pixels. These regions should be representative 
and can be made up of different not connected parts. The 
unsupervised algorithm is an improvement of the supervised 
ones and is completely automatic. 

To estimate the parameters of the class conditional densities, 
a new expectation maximization (EM) algorithm was devel- 
oped. Details are described in appendix. When segmenting 
small sub-scenes of the image, a simplified data model with 
only one Gamma function per class can be used. 

The prior used to impose local homogeneity is a multi- 
level logistic (MLL) model Markov random field (MRF) [?], 
with 2nd order neighborhood. To infer the prior smoothness 
parameter controling the degree of scene homogeneity, we 
develop an EM algorithm that uses loopy believe propagation 
(LBP) [?]. We have also exploited different classic estimation 



methods, namely the least squares fit (LSF) and the coding 
method (CD) (see [?] for details) for comparison. 

To infer the labels, we adopt the maximum a posterior 
(MAP) criterion, which we implement efficiently and exactly 
with graph-cut techniques [?]. 

Although the segmentation of oceanic dark patches is 
typically based on the assumption of two classes, we have 
generalized the problem to an optional number of classes. The 
underlying integer optimization problem is now attacked with 
the graph-cut based a-Expansion algorithm [?]. 

To evaluate the accuracy of the proposed algorithms, differ- 
ent simulations addressing both the referred Gamma model as 
well as intensity images corrupted with Gaussian noise have 
been carried out for error rate assessment. 

The algorithms have also been applied to real SAR images 
containing well documented oil spills. For doing so, the scenes 
are divided in tiles and segmented individually. 

A. Related Work 

Related approaches to the problem of oil spill segmentation 
are built on off-the-shelf segmentation algorithms such as the 
adaptive image thresholding and the hysteresis thresholding. 
Entropy methods based on the maximum descriptive length 
(MDL) and wavelet based approaches have also been pro- 
posed. Another recently proposed segmentation methodology 
applies Hidden Markov Chains (HMC) to a multiscale repre- 
sentation of the original image. Hereby the wavelet coefficients 
are statistical characterized by the Pearson system and by the 
the generalized Gaussian family [?]. When other SAR products 
are available, for example polarimetric data, other methods 
have been described in the literature, like constant false alarm 
rate filters [?]. 

An example of an elaborated adaptive thresholding tech- 
nique is provided by [?]. In this method, an image pyramid 
is created by averaging pixels in the original image. From 
the original image, the next level in the pyramid is created 
with half the pixel size of the original image. A threshold 
is then computed for each level based on local estimates of 
the roughness of the surrounding sea and on a look-up table 
containing experimental values obtained from a training data 
set. 

Hysteresis thresholding has been used as the base for 
detecting oil slicks in [?]. The method includes two steps: 
applying a so-called directional hysteresis thresholding (DHT) 
and performing the fusion of the DHT responses using a 
Bayesian operator. The MDL technique, which basically con- 
sists in applying information theory in order to find the 
image description which has the lowest complexity, has been 
applied in [?] to segment speckled SAR images, namely those 
containing oil slicks. This segmentation method describes the 
image as a polygonal grid and determines the number of 
regions and the location of the nodes that delimit the regions. 
The two-dimensional wavelet transform, used as a bandpass 
filter to separate processes at different scales, has also been 
adopted to oil slick detection in the framework of an algorithm 
for automated detection and tracking of mesoscale features 
from satellite imagery. [?]. 



B. Contributions 

We approach oil spill segmentation using a Bayesian frame- 
work and a multi-level Logistic (MLL) prior. Several methods 
in the same vein have been proposed since the seminal work 
of Geman and Geman [?], see e.g., [?]. Applications of these 
ideas in the segmentation of SAR images can be found, e.g., 
in [?], [?], [?]. 

The main contributions of this work to the state-of-art in oil 
spill segmentation are the following: 

• the development of an EM algorithm to estimate the 
parameters of a mixture of a pre-defined number of 
Gamma distributions, in order to model the intensities 
in a SAR image 

• the development of an EM algorithm using LBP to 
estimate the smoothness parameter in the MRF used as 
pior in our framework 

• the application of recent graph-cut techniques for solving 
the energy minimization problem that arises from the 
followed Bayesian methodology. 

• the design of a semisupervised algorithm for oil spill 
segmentation supported on the tools referred to above. 

C. Paper Organization 

The article is organized as follows: Section 1 introduces 
the problem, with references to related work and giving the 
main contributions of the present work; Section 2 overviews 
the Bayesian methodology that builds the base to the proposed 
algorithms. In addition, it briefly reviews the concept of the a- 
Expansion technique that has been implemented to generalize 
the methodology to an optional number of classes; Section 3 
describes the used parameter estimation techniques; Section 4 
describes the supervised segmentation algorithms; Section 5 
describes the unsupervided algorithm and Section 6 presents 
results of segmenting simulated and real images applying the 
algorithms proposed in Section 4 and 5; Finally Section 7 
contains the main conclusions and future work remarks. The 
article also includes an appendix where the referred to EM 
algorithm, developed to estimate the parameters of the class 
conditional densities of the Gamma mixture data model, is 
described. 

II. Problem Formulation 

A. Bayesian Approach 

let C := {l,...,c} be a set of c classes and V := 
{1,2, . . . , A/"} be the set of N pixels (sites) where mea- 
surements y := {^1,^2, . . . ,^Ar}, the SAR intensities, are 
available. A labeling x := {xi, X2, . . . xn} is a mapping from 
V to £, i.e., it assigns to each pixel p e V 3. label Xp G C. 
Any labeling x can be uniquely represented by a partition of 
image pixels P = {Pi\l G jC}, where Pi = {p e V\xp = 1} is 
the subset of pixels to which the label / have been assigned. 
Since there is an one-to-one correspondence between labeling s 
X and partitions P, we use these notions interchangeably. By 
applying a segmentation algorithm to the image y, we get 
X := {xi, ^2 . . . Xn}, where Xi is the inferred label for pixel 
i e V. 



B. Observation Model 

In our problem formulation, we assume conditional inde- 
pendence of the measurements given the labels, i.e., 
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p{y\x) = Y[p{yi\xi) = n n ^(^^1^^)' 
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where p(-|0O is the density corresponding to class / and (j)i 
the correspondent vector of parameters. The adopted density 
is a finite Gamma mixture given by 
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p{y^\^') = Y.<P(y^\^'s)^ 
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where K is the number of Gamma modes in the mixture, 
i indexes the pixel, and, for the class /, O^ is the vector of 
parameters of the Gamma mode s, and a[ is the a priori 
probability of mode s. We denote 6^ := {0[, . . . , 6>^), a^ := 
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Given that p{yi\0\) is Gamma distributed, we have then 
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y°' 'exp(-Aiyi), t/i > 0, (3) 



where ^^ := (a^,A^). The mean and variance of a random 
variable with the density dSji is, respectively, a^/Ag and 

The mixture parameters 0^ are estimated from the data 
by applying the EM Gamma mixture estimation algorithm 
described in Appendix. The procedure is the same both for 
the supervised and for the unsupervised algorithms. 

We now make a brief comment on our choice of the Gamma 
mixture for modeling the observation densities. Under the 
assumption of fully developed speckle, the complex radar 
amplitude is zero-mean circular Gaussian distributed [?]. The 
average intensity computed over a number of independent 
random variables with the same density is, thus. Gamma 
distributed. It happens, however, that in sea SAR imaging 
one or more of the above assumptions may fail, rendering 
the Gamma density a poor model for SAR intensity [?]. A 
line of attack to obtain better models is to use more flexible 
parametric families, such as the Pearson System [?, Ch. 4.1] 
(see also [?]). However, in our problem each class density 
exhibits much more variability than that of accommodated 
by the the Pearson System. We have very often, for exam- 
ple, multi-modal densities. We should resort, therefore, to 
a mixture of Pearson System densities, what would lead to 
complex learning procedures. We have experimentally ob- 
served, however, that the Gamma mixture yields very good 
fittings for real SAR histograms, obtained with a moderately 
complex learning algorithm. For this reason, we have adopted 
the Gamma mixture model. 

C. Prior 

A second assumption we are making is of local smoothness 
of the labels in a statistical sense. It is more likely to have 
neighboring sites with the same label than the other way 
around. We model this local smoothness with a second order 



MRF, P {x)\ more specifically, by an MLL model (Ising model 
in the case of two classes). The Markov property assumes that 



p{xi\ XjJ e V) =p{xi\ XjJ eMi), 



(4) 



where Mi is the set of neighbors of site i. If 
p{xi\ Xj.,j ^ J\fi) > 0, then the Hammersley-Clifford 
Theorem states that p{x) has the form of the Gibbs 
distribution 

p(x) = iexp-^(-), (5) 

where Z is the so called partition function and U is the energy 
function 



U{x)= ^Vci{a 



(6) 
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where C is the set of cliques and Vd (x) is the clique potential 
defined over clique d. In this work, we have pair- wise cliques 
defined on a second-order neighborhood (8 pixels). That is, 
C = {{ij) : i e Mj, j eMi, i> j). 

In these conditions, the MLL clique potentials, in the 
isotropic case, is given by 



Vcl {Xr, Xs) = -PS {Xr - Xg) , 



(7) 



where r, 5 G cl, S{x) := Io{x) is the indicator function of 
set {0}, and parameter (3 > controls the degree of scene 
homogeneity. 



D. Maximum a Posteriori Estimate 

The posterior of the labeling given the observed data is 

p{x\y) (X p {y\x) p (x) . (8) 

In order to infer x, we adopt the MAP criterion. This amounts 
to maximize the posterior density of the labeling given the 
observed data: 



X = argmaxp {x\y) , 



(9) 



and is equivalent to minimize the negative logarithm of ([8]). In 
this sense, we may rewrite the problem in the following way: 



where 



X = argmin£^(xi, . . . ,Xn), (10) 



E{xu...,xn) = -logp(x|y)+c^^ (11) 



and c^^ denotes an irrelevant constant. From equation ([8]), we 
have 
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E(xi,...,x^) = ^E^(x,)+ ^ E'^^{x,,xj), (12) 



with 



E'{x,) = -\ogp{y,\x,) (13) 

E'^\xi,Xj) = -f36{xi-Xj). (14) 



E. Energy Minimization 

As already stated, we are concerned with the minimization 
of E {xi^. . .xn) given by ([T2b . which we term energy. 
For two classes, {i.e., c = 2), the global minimum of 
E {xi^. . . xn) can be computed exactly by applying the graph- 
cut algorithm described in [?]. This is a consequence of energy 
being graph-representable, i.e., £^*'-^(0,0) + £^*'-^ (1, 1) < 
E^^^ (0, 1) + E^'^ (1, 0) (for details see [?]). For more than two 
classes, the solution of Q can be approximately computed 
by the a-Expansion technique [?], also based on graph-cut 
concepts. This algorithm finds the local minimum of the 
energy within a known factor of the global minimum. 

We now give a brief description of the a-Expansion al- 
gorithm. Given a label a, a move from a partition P (with 
correspondent labeling x) to a new partition P' (with corre- 
spondent labeling x') is called an a-Expansion if Poc C P'^ and 
P/ C Pi for any label / ^ a. In other words, an a-Expansion 
move allows any set of image pixels to change their label to 
a. The algorithm cycles through the labels in V in some fixed 
or random order and finds the lowest a-Expansion move from 
the current labeling. If the expansion move has lower energy 
than the current labeling, then it becomes the current labeling. 
The algorithm terminates when no a-Expansion move exist 
corresponding to a local minimum of the energy. 

III. Estimation of Parameter p 

In this work, we consider three different techniques to 
determine the smoothness parameter /3: a new EM method 
hereafter introduced that uses loopy belief propagation (LBP) 
and the classical LSF and CD methods. Because LSF and CD 
assume the existence of labeled data, we have conceived an 
iterative labeling-estimation scheme, which alternates between 
a labeling step and an estimation step until the convergence of 
P is attained. On the contrary, the EM estimation algorithm, 
that we have called "Loopy-/3-Estimation", is a one-shot 
technique. In the following sections, we briefly review the 
LSF and CD methods and provide a detailed description of 
the Loopy-/3-Estimation method. In this section, the class 



parameters 



, . . . , ^^) are assumed known. 



A. Least Squares Fit 

This procedure for parameter estimation in MRF, described 
in detail in [?], is based in the following equation that holds 
for the MLL model, for every pixel p, with neighborhood A/^, 
and for every label pair Xp^Xp e C : 



/3 [n{xp) -n{xp)] 



log 



'p{xp\xxp,y) 



^p{x'p\xM^,y) ) 

- {E^ (x;) - EP {xp)) , (15) 

where n{xp) is the number of pixels in the set JVp with 
the same label as Xp and xj^f^ := {xi^ i G JVp}. We use 
histograms to estimate the joint probabilities p {xp\xj\fp^y) 
and p (^Xp\x^p^y) in ([T5]) : assuming that there are a total of 
M distinct 3x3 blocks in the image lattice with a given label 
configuration xj^f^, then we take 

y^ {xp\xj^p,y) 



where H {xp\xXp^y) is the number of times that a particular 
3x3 configuration {xp\xXpjy) occurs. The expression is then 
evaluated for a number of distinct combinations of Xp , x^ and 
xj\fp in order to obtain an over-determined linear system of 
equations that is solved in order to /3. 

B. Coding Method 

In this method the key idea (see [?]) is to partition the set 
V into sets V^^\ called codings, such that no two sites in one 
-pik) ^^Q. neighbors. In the present work, the neighborhood ATp 
is of 2nd order thus yielding four codings. As the pixels in 
-pik) ^^Q. j^Q^ neighbors, the variables associated with these 
pixels, given the labels at all other pixels, are mutually 
independent under the Markovian assumption. The following 
simple product is thus obtained for the likelihood: 

p(^)(x|/3,y)= \[ p{xp\xMp,P,y), (17) 

pev^^^ 

Maximizing (fTTl) in order to P gives the coding estimate 
P^^\ Although it is not clear how to combine the results 
optimally, the arithmetic average, as suggested in [?], is an 
intuitive scheme that was adopted in this work. 

C. Loopy- j5 -Estimation 

We are seeking 13ml = argmax/3p(?/|/3), the ML estimate 
of the smoothness parameter /3. Based on the fact that the 
marginal density p{y\P), the so-called evidence, is a sum over 
the missing labels x, i.e., 



p{y\P) = ^p{y^x\P) 

X 

= ^p{y\x)p{x[ 



(18) 



we develop an EM algorithm [?] for the ML estimation 
of the parameter /3. The EM algorithm alternates between 
two steps: the E-step computes the conditional expectation 
of the logarithm of the complete a posteriori probability 
function, with respect to the missing variables, based on the 
actual parameter value; the M-step updates the value of the 
parameter, by maximizing the expression obtained in the 
E-step with respect to that parameter. We now derive the 
E-step and the M-step. 



p[xp\xM^,y) = 



M 



(16) 



E-step: 










Q{l3;l3t) = E[logp{y,x\l3)\y,l3t] (19) 


= E[logp{y\x)\y,/3t] (20) 


+E[\ogp{xmy,pt]. (21) 


Recalling that the MLL prior is given by 


pN/3)=2(^)exp 
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with 
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Fig. [T] depicts the graph that represents our pairwise MRF 
for computing the two-node behefs. In the square lattice, 
ipij {xi , Xj ) stands for the interaction potential that penalizes 
every dissimilar pair of neighboring labels, and (j){xi^yi) = 
p{yi\ Xi) represents the statistical dependency between the 
labels Xi and the measurements yi. For computing hij (xi^Xj), 
(j){xi,yi), is set to a constant value, making the result inde- 
pendent of the data values y. 

We solve ([28l) by a line search type algorithm, ensuring that 



< 0, thus corresponding to a maximum of Q. 



Fig. 1. Lattice representing the pairwise MRF. 



we obtain, up to an irrelevant constant. 



Q(/3;A) = -log^(/3) + /3 Y. E\b{x,-x,)\y,^t\. 



M-step: 



(3t+i = arg max Q (/?, A) . 



(24) 
(25) 



The stationary points of Q are the solution of 

dQ _ d\og[Z{/3)] 
dp dp 



(26) 



+ X! X! Hxi-Xj)p{xi,Xj\ y,/3t))=0. 



By introducing expression (1231) into (1261) and if we consider 
that 5 {xi — Xj) takes non-zero values only for equal labels, 
we obtain 



dQ 

dp 



i,jGcl k—1 

-p{xi = k,Xj = k\ Pt) = 0. 



(27) 



Since computing exact marginal distributions is infeasible 
in our case, we replace them by pseudo-marginals using BR 
This approach has been successfully applied in problems of 
approximate parameter learning in discriminative fields [?]. 
BP is an efficient iterative algorithm in which local messages 
are passed in graphical models. For singly-connected (loop 
free) pairwise MRFs, the two-node beliefs will correspond 
to the exact two-node marginal probabilities. In our case, 
however, the graph that corresponds to the MRF contains 
cycles, preventing the basic BP algorithm to be applied. 
We resort to a slightly modified version called loopy belief 
propagation. In practice, this algorithm has often delivered 
good results [?]. We approximate the marginal probabilities 
with the two-node beliefs bij (xi^Xj) and b^j {xi^Xj). These 
will provide approximations respectively for the marginals 
p{xi^Xj\l3t) and p{xi^Xj\y^ Pt). By doing so, the M-Step of 
our EM algorithm is given by the sum of the differences 
between the two node beliefs that take the evidence y into 
consideration and those that do not make use of y\ 



^ i,jed fe=l 



■ bij {k, k) = 0. 



(28) 



D. A few Remarks about the Vector of Parametres (j) 

In the previous section, we have considered that the class pa- 
rameter vector (j) is known and only the smoothness parameter 
P is to be inferred. However, still using the beliefs computed 
by LBP, the vector (j) could have been inferred simultaneously 
with p by including in the Q function the additive term (l2Qb 
corresponding to the class densities. We would have obtained 

Q{^',^uPt) = E[\ogp{y\x)\y,(l)upt] 

N c 

= X]X]logp(?/i|^^)p(xi =l\y,(t)t,pt), 

where the probabilities p{xi = l\y^(j)t^Pt) are given by the 
LPB method. The M-step would consist then in two decoupled 
maximizations; one with respect to /3 and another with respect 
to (j). This approach is, however, beyond the scope of this 
paper. Nevertheless, we present in Section |Vl an unsupervised 
algorithm, which is suboptimal but faster than the LBP based 
EM algorithm aimed at the inference of both the parameters 
P and (j). 

An alternative to the proposed EM scheme based on LPB 
is using Monte Carlo techniques to cope with the difficulty 
in computing the E-step and the partition function Z [?], [?], 
[?]. Supported on the performance of the LBP, we believe, 
however, that the EM scheme based on LPB is, for the 
present problem, much more faster than the Monte Carlo based 
techniques. 

IV. Supervised Segmentation 

In this section, we introduce two supervised algorithms 
aimed at the segmentation of sea SAR images. In both algo- 
rithms, the first step is the estimation of the class parameters 
used in the data model. This is done by asking the user to de- 
fine representative regions of interest (ROI) in the image. Once 
the regions are defined, different approaches may be followed: 
if a Gamma mixture is assumed for the observed data, as in 
the case of SAR images, the EM class parameters estimation 
algorithm described in Appendix is applied to infer the class 
conditional densities. If we are segmenting small sub- scenes 
and the radar range spreading loss has been compensated in the 
underlying SAR images, one single Gamma function per class 
often provides a good modeling for the SAR intensities. In this 
case, a common ML Gamma estimator is used instead of the 
EM procedure. After this step, the data model is considered 
to be known and is used thereafter. 



The pseudo-code for two supervised algorithms is presented 
in Algorithm 1 and Algorithm 2. Algorithm 1 is of generalized 
likelihood type [?] implementing an iterative labeling scheme 
with two steps being performed alternately: the /3 -Estimation 
and the segmentation. Algorithm 2 is a one-shot procedure that 
performs /3-Estimation using the Loopy /3-Estimation method 
described in Section IIII-CI 

Algorithm 1 Supervised Segmentation Using LSF/CD (3- 

Estimation 

Require: Initial parameter (3 = Pq and estimated class pa- 
rameters (j) 
Compute eLabel = E^{xp) for every pixel p using ^ 



while 



A/3 



< 5 OX nrlterations < NrlterationsMax do 



Compute X 
Compute P 

end while 

return (x,/3 



a-Expansion(;^, eLabel) 
l3-Estim3.tion{x^ eLabel) {LSF or CD}. 



Algorithm 2 Supervised Segmentation Using Loopy-/3- 

Estimation 

Require: estimated class parameters cj) 
1: Compute eLabel = E^{xp) for every pixel p using (p 
2: Compute /3 = Loopy-/3-Estimation(eLa6eO 
3: Compute X = GraphCut-Segmentation(/3, eLabel) 
4: return (x,^) 



V. Unsupervised Segmentation 

The unsupervised method is an improvement of the super- 
vised algorithms described in the previous sections. The main 
difference is that the data model is not considered to be known 
but is also iteratively estimated along with the smoothness 
parameter and the segmentation. The scheme needs a rough 
initialization of the data model parameters. We computed this 
initialization by fitting an EM Gamma mixture to the complete 
data. 

As our EM algorithm automatically eliminates unnecessary 
modes, we start with an overestimate of K, the number of 
modes in the mixtures. In our experiments with real data, the 
maximum number of modes we got was four. 

Then, different strategies are possible. In the experiments 
reported in the next section, when considering c = 2 classes, 
we have assigned the mode with a lower mean value to one 
of the classes and the remaining modes to the other class. 

In each iteration, we compute the ML estimate of the 
vector (j) based on the previous segmentation, compute the new 
segmentation, and finally the new value of f3. The algorithm 
may use any of the three /3— parameter estimation methods and 
is applicable to an optional number of classes. The algorithm 
stops when both the parameter (3 and the class parameters 
converge to stable values. 



Algorithm 3 Unsupervised Segmentation 

Require: arbitrary parameter /3 = /3o, initial class parameters 
(j) = (j)Q (EM Gamma Mixture Estimation) 
1: Compute eLabelo = EP{xp) for every pixel p, using 0o 
2: Compute initial labeling 

X = xo = a-Expansion(/3o, eLabelo) 
3: for Stop criterium is not met do 
4: Compute (/) =ML-Estimation(f) 
5: Compute eLabel = E^{xp) for every pixel p, using ^ 
6: Compute X = Qf-Expansion(/3, eLabel) 
7: Compute $ = /3-Estimation(x, eLabel) {use LSF, CD 

or Loopy}. 
8: end for 
9: return (x,/3,^j 



VI. Results 

This section presents results of applying the proposed 
methodology to simulated and real SAR images, as well 
as some considerations regarding time complexity of the 
algorithms. 

A. Simulations 

We have performed simulations corresponding to Gamma 
data terms with two classes, and evaluated the overall accura- 
cies of the obtained segmentations using Algorithms 1 , 2 and 
3. We also illustrate in detail the EM Gamma mixture esti- 
mation algorithm, by applying it to an example of simulated 
data. 

1) Segmentation Results with Gamma Data Term: In order 
to compare the performance of the three proposed methods, 
we have tested Algorithms 1, 2 and 3 on simulated images 
generated by adding Gamma noise to ground-truth images 
containing two classes. We have used the ground-truth de- 
picted in Fig.|2](a) and added noise with Gamma distributions 
having mean values of five and nine. The parameters of the 
distributions were choosen in order to obtain increasing values 
of variance a'^, corresponding to noisier images. We have then 
applied Algorithm l,with LSF and CD methods. Algorithm 2, 
and Algorithm 3 with Loopy-/3-Estimation. Furthermore, for 
comparison, the images were also segmented tuning the beta 
value manually and for the case that no prior is used (13 = 0). 
Fig. [2] shows the segmentations obtained for a = 2.6. As we 
can see from /3 = this is a hard problem, on which LSF and 
CD fails the inner structures, but Loopy-/3-Estimation provides 
better results, both in the supervised as in the unsupervised 
method. This behavior is further confirmed by the rank in Fig. 
[51 which gives the OA obtained for images with different a 
values (corresponding to images with more noise) for the six 
different segmentation processes referred above. 

We have also tested the proposed algorithms in images 
generated by adding Gamma noise to a ground-truth similar 
to an oil patch (see Fig. \^. Two images. Image A and Image 

B, with the histograms depicted in Fig. IH corresponding to 
different segmentation difficulty levels are segmented. The best 
segmentation possible in this framework, achieved by tuning 
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Fig. 2. Ground-truth and results of different segmentation processes for an 
image with Gamma noise having mean values of five and nine and with a a 
value = 2.6. Notice the good performance of Algorithm 2 and Algorithm 3, 
implementing the Loopy-/3-Estimation. 



the /3 value manually, and the segmentation obtained with 
no prior information (setting /3 = 0) are also displayed for 
comparison. Fig. [5] shows the results obtained for Image A 
and Fig. [6] the results obtained for Image B. For Image A, the 
segmented images using Algorithm 1 , with CD and with LSF, 
are not shown, as they are almost equal to the image segmented 
with Algorithm 2. For Image B, Algorithm 3 has not provided 
good results and the segmentation is not displayed. The bad 
performance of Algorithm 3 in this case arises from the not so 
good initialisation of the class parameters, due to the complete 
overlapping from the oil and the water histogramms. 

2) EM Algorithm for Gamma Mixture: In this subsection 
we illustrate the behavior of the EM algorithm designed to 
infer Gamma mixtures. For details on the theoretical issues, 
we refer to the Appendix. In Fig. [71 we can see the ground- 




Fig. 3. OA for images corrupted with Gamma noise with increasing a values: 
TM = for /3 tuned manually; LE = Algorithm 2, using Loopy Estimation; 
LSF = Algorithm 1 using LSF; CD = Algorithm 1 using CD, NP = no prior 
and UNS = Algorithm 3 using Loopy Estimation. 
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Fig. 4. Gamma Model: (a) Ground-truth, (b) Histogram of Image A, 
(c) Histogram of Image B, (d) Image A, (e) Image B with superimposed 
delimiting line, for better visualization 
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(e) 

Fig. 5. Segmentation results for image A: (a) No prior information, (b) f3 
value tuned manually, (c) Segmentation using Algorithm 3 with Loopy-^- 
Estimation, (d) Segmentation using Algorithm 2 (Algorithm 1 provided the 
same results with only 0.1% difference in OA), (e) graph depicting the ^ 
estimation in Algorithm 2. 



truth used for simulating the image and the generated image, 
according to the densities shown in Fig.[8l The two classes, oil 
and water have been modeled respectively by a mixture of two 
Gamma functions and a mixture of three Gamma functions. 
For this particular case, we have selected a ROI containing 
178 pixels for representing the water and a ROI containing 142 
pixels for representing the oil. After 20 iterations of the EM 
algorithm, we obtained the approximations for the probability 
distributions depicted in Fig. [9] and Fig. \TU\ The obtained 
results are quite reasonable for the small sample sizes used. 

B. Real Images 

We have tested the proposed methodology with three real 
SAR images for the special purpose of oil spill detection. For 
this type of application, an unsupervised algorithm is more 
indicated, and so we decided to apply Algorithm 3. Because 
Loopy Estimation has proven to be an effective method in 
the simulations, we have choosen Algorithm 3 with Loopy 
Estimation. Nevertheless, in order to have a comparision 
unsupervised approach versus supervised approach, we also 
applied Algorithm 2 to two of the images. In oil spill detection. 
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Fig. 6. Segmentation results for Image B: (a) No prior information, (b) (3 
value tuned manually, (c) graph depicting (3 estimation using Algorithm 2, (d) 
Segmentation using Algorithm 2, (e) Segmentation using Algorithm 1 with 
LSF and (f) Segmentation using Algorithm 1 with CD. 





(a) 



Fig. 7. (a) Ground-truth used for simulating an oil spill. Black represents 
oil and white water; (b)Simulated SAR image 



the number of classes is typically set to c = 2, although more 
classes can be considered if we are interested in distinguishing 
other phenomena occuring at the same time in the region of 
interest. The application of the algorithms is straightforward: 
the key idea, like in most state-of-the-art oil spill detection 
methods (see for example [?]), is to partition the image in 
tiles and run the algorithm separately for each part. This step 
is preceeded by the application of a landmask to the image, 
what can be done using external coastline information or by 
adopting some coastline self-extraction procedure. After the 
segmentation of each tile, a procedure for grouping patches 
detected on the tile borders should be carried on. Another 
possibility to increase segmentation coherency in the borders 
is to define overlapping tiles or to force continuity to some 
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Fig. 8. Probability functions used to generate the simulated image with 
superimposed histogram of generated data set. A three modes function for 
water and two modes function for oil was used. 
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Fig. 9. True and estimated class densities for oil. 
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Fig. 10. True and estimated class densities for water. 




Fig. 11. ERS-1 image from the Sicily Channel, Italy, acquired on 30th 
January of 1992. The smaller and larger squares are sub-scenes that have been 
used respectively to estimate the wind direction and to apply our algorithm. 




m 



(b) 



Fig. 12. Closer look to the sub-scenes of the ERS-1 image. Part (a) is used 
for wind direction estimation; the estimated direction has been overlapped as a 
green arrow. Part (b) contains an oil spill. The coloured regions correspond to 
the ROI's selected for class parameter estimation in the supervised algorithm 



degree on the estimated class and/or smoothness parameters 
from one tile to the next. Nevertheless, at the moment we 
are not doing this and these are considered future possible 
improvements to our methodology. We segmented oil spills 
contained in three different scenes, described in the following 
subsections. 

1) Segmentation of an ERS-1 sub-scene: We have seg- 
mented part of an ERS-1 image from the Sicily Channel, 
Italy, that has been acquired on the 30* January of 1992. 
This image is referred in the ESA web pages regarding oil 
slicks http : //earth . esa . int /ew/oil_slicks/ and 
contains three oil slicks, along with information regarding 
wind direction and intensity and existence of ships, ship wakes, 
natural oil films and currents. 

Fig. [TT] provides a quicklook of the scene with two 
squares superimposed: the larger representing the part to be 
segmented and the smaller representing a part used for wind 
estiamtion. Fig. [121 provides zooms of the referred squares. 
By applying the Radon Transform to the smaller square, the 
estimated direction has been calculated and is depicted in the 
image. The direction was consistent with the measured value 



10 



4.5 

4 
3.5 

3 
2.5 

2 
1.5 

1 
0.5 






0.4 0.5 0.6 0.7 0.8 0.9 



Fig. 13. Fitting of a mixture of Gammas to the data from the red square in 
the ERS-1 image. 




Fig. 14. Segmentation of ERS-1 subscene containing oil with (3 
estimated using Algorithm-3 with the Loopy-^ -Estimation method. 



1.4, 



reported in the url site and is at the origin of the well-known 
"feathering" effect that can be observed in this linear spill. 
We have computed the backscattering values of the image, by 
performing calibration using the ESA provided BEST software 
Chttp : //earth. esa. int /services /best/) and 
then applied Algorithm-3 with Loopy-/3-Estimation. The 
result of the Gamma mixture estimation, in the initialisation 
of the algorithm, is depicted in Fig. [131 After only three 
iterations, both the class parameters and the smoothness 
parameter have converged. A /3 value equal to 1.4 was 
estimated. The segmentation is displayed in Fig. [141 We also 
applied Algorithm-2: we selected two ROI's (160 pixels for 
water and 77 for oil) in the image (shown in Fig. [12]) and 
computed ML Gamma Estimators for the two classes (see 
Fig. [15]). We then applied the algorithm and obtained an 
estimated P value equal to 1.44. The segmentation result is 
given in Fig. [16] and is similar to the result obtained with 
Algorithm-3. 

For comparison, we also provide the segmentation obtained 
with no prior (corresponding to /3 = 0) in Fig. [TT] In practice, 
the estimated value has proved to deliver a good segmentation. 
When lower /3 values were used, the result was noisier and 
for higher /3 values, the details of the spill disappeared. 
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Fig. 15. Class Parameters Estimation for Algorithm-2. 




Fig. 16. Segmentation of ERS-1 subscene containing oil with ^ 
estimated using Algorithm-2. 
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Fig. 17. Segmentation of ERS-1 subscene containing oil, with /3 = 0. 




Fig. 18. (a) ASAR image fragment with coloured regions for class parameter 
estimation in Algorithm-2; (b) Segmentation result applying Algorithm-2; (c) 
Segmentation result applying Algorithm-3 with Loopy-/3-Estimation. 




Fig. 19. Display of the ASA WSM full resolution image of the Prestige 
accident occured in November 2002 in Galicia. 



2) Segmentation of an Envisat ASAR IM sub-scene: 
We have also applied Algorithm 3 to a fragment of an 
ASAR Image Mode image, acquired on 19 July 2004, in 
the ocean between Cyprus and Lebanon. The fragment 
contained an occurred oil spill of circa 10 km's (centered 
on ^ 33°N,33°E39') that was documented on the EC 
Oceanides project (a project in the framework of the Europe 
's Global Monitoring for Environment and Security initiative) 
database. After six iterations, we achieved convergence of 
the P and values. The segmentation result is given in Fig. 
[TSl corresponding to a estimated /3 = 1.83. The unsupervised 
segmentation result is compared with the one provided by 
Algorithm-2, where the user provided ROI's for water and 
for oil that were used to estimate the class parameters at 
the beginning of the process (see Fig. [18]). In this case a 
/3 = 1.75 was obtained. 

3) Segmentation of an Envisat ASAR WSM sub-scene: In 
order to demonstrate the viability of applying Algorithm 3 to a 
whole ASAR_WSM scene, we have run it over the very well 
known image of a confirmed oil spill, namely the Prestige 
case (see Fig. [19]). This accident took place in November 
2002, in Galicia (Spain), when a tanker carrying more than 
20 million gallons (around 67,000 tons) of oil split in half 
off the northwest coast of Spain on 19 November 2002, 
threatening one of the worst environmental disasters in history. 
For segmenting the image, we have first partitioned it in tiles of 
600x600 pixels each and then applied Algorithm-3 to each tile 
independently. No post-processing of the borders or "clean-up" 
operations (like for example morphological operations) have 
been carried out. Fig. [20] shows the segmented image, made 
up of the concatenation of the segmentation of the individual 
tiles. As we can see, the result can be considered in general 
very good, with only some discrepancies located on the tiles' 
borders. For initialising the algorithm, the EM Gamma mixture 



estimation procedure was applied with four modes {K = 4). 
In fact this number is enough to provide a good fitting of the 
intensity data of the SAR image. On the other hand, our EM 
algorithm allows us to start with a higher mode number, and 
decreases this number automatically. We depict the results of 
this fitting on one of the tiles, shown in Figure [23] by showing 
the histogram of the tile with the superimposed estimated 
Gamma mixture (see Fig. [21] (a)). As we have explained, the 
mode corresponding to the lower mean value is assigned to 
the oil class and the others to the water class, providing an 
initialisation to the class parameters. The estimations of these 
parameters are actualised along the algorithm and, after nine 
iterations, we obtain the distributions shown in Fig. [21] (b). 

To fully demonstrate the possibilities of our algorithm, 
we again run it on the tile shown in Figure [23] but this 
time setting the number of classes to c = 3. By doing so, 
we hope to be able to segment a third ambiguous zone, 
corresponding to intermediate radiometry levels and probably 
due to atmospheric conditions originating a front. In this case 
we choose to fitt a Gamma mixture of c modes to the data 
for initialising the class parameters, as depicted in Fig. [22] 
With this initialisation, the first obtained segmentation using 
Algorithm 3 is displayed in Fig. [24] after only nine iterations, 
all parameters have already converged to the segmentation 
given in Fig. [25] When comparing this segmentation with 
the one provided by a state-of-the-art algorithm, namely by 
multiscale HMC model in [?], we consider to have obtained 
a very good result. 

4) Time considerations: The proposed algorithms display, 
in practical cases, a computational complexity of 0{N), with 
N being the number of pixels in the image. In fact, the most 
time-consuming steps in the algorithms are the a-Expansion 
and Graph-Cut segmentation routines. These use the max- 
flow code implementation referred in [?] that has complexity 
0{n'^-^) in the worst-case, however, in most situations is 
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Fig. 20. Display of the segmentation results of applying Algorithm 3 with 
Loopy-/3-Estimation to the image displayed in Fig. [19] 





Fig. 22. Histogramm of the data values (normalised) with superimposed 
estimated Gamma mixture when the number of classes is set to c = 3 




Fig. 23. Display of one of the individual tiles in which the image in Figure 
[191 has been divided for segmentation. 
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Fig. 21. Results of applying the EM Gamma mixture estimation procedure to 
the tile depicted in Figure [23] (a) Histogramm of the data values (normalised) 
with superimposed estimated Gamma mixture; (b) Distributions corresponding 
to the oil class (left) and water class (right) after nine iterations of Algorithm 
3: Histogramm and superimposed fitting. 




Fig. 24. Display of the initial segmentation result of applying Algorithm 3 
with Loopy-^-Estimation and three classes to the tile depicted in Figure [23] 




Fig. 25. Display of the final segmentation result, after nine iterations, 
of applying Algorithm 3 with Loopy-/3-Estimation and three classes to the 
bottom left tile of Fig. 17. 



0{N). Both Algorithm 1 and Algorithm 3 converge in few 
iterations and, as an example, for images of the dimension used 
in the first simulation (64x64) (see lVI-ATTl) , only a few seconds 
are needed to run the algorithms. For a 600x600 pixel tile, like 
those in the Prestige example in lVI-B.3l Algorithm 3 is taking 
^ 1 minute, but this time performance could be improved by 
a more efficient implementation of our procedure, namely by 
improving the implementation of the Loopy-/3-Estimation and 
by migrating our matlab code fully to C-code. 

VII. Conclusions 

The results of applying the proposed methodology to simu- 
lated images with Gamma data models and to real SAR data 
are promising. The developed EM Gamma mixture estimation 
algorithm, when incorporated into the proposed segmentation 
algorithms, has proved to be an efficient tool for data modeling 
in SAR intensity images. 

With the supervised algorithms, high OA accuracies have 
been achieved even for simulated images with high levels 
of noise. In general the Bayesian approach has resulted in 
an OA increase in the segmentation process between 10 and 
15%, when compared to the segmentation using no prior 
information. Hereby, Algorithm 1 with LSF and Algorithm 1 
with CD methods provided similar results, with a performance 
close to that obtained by setting the /3 value manually. On the 
other hand, the LSF estimation procedure seems to become 
less reliable for noisier images. When compared to Algorithm 
2, using Loopy-/3-Estimation, we see that this last method 
allways provides equal or better results, outcoming Algorithm 
1 . Another advantage of Algorithm 2 is that it is usually faster 
than Algorithm 1, being a one-shot process. Both algorithms, 
by introducing prior information into the segmentation pro- 
cess, increase the OA significantly. As a conclusion. Algorithm 
2 should be preferred to Algorithm 1 , when using a supervised 
method for oceanic sar images segmentation. 

Algorithm 3, totally unsupervised, has been conceived as 
an improvement to Algorithm 1 and 2. When testing it on 
simulated data the obtained results were very good, although 



slightly worst than the supervised ones, as expected. In gen- 
eral, after a few iterations, the class and smoothness param- 
eters converge to stable and meaningful values. By applying 
Algorithm 3 to real images containing documented oil spills, 
the segmentation has been considered successful! We could 
segment both linear and patch type oil spills. Furthermore, the 
applicability of the method to segment whole scenes, as well as 
to segment more than two classes, has also been demonstrated 
in a well known image from the Prestige accident. As a 
conclusion, we believe the presented methods are suitable to 
be used for segmenting oceanic SAR images. In particular for 
oil spill detection. Algorithm 3 seems to be a suitable method. 
Based on observations that oil spills in the ocean are often 
dragged by the wind and align more or less perpendicular to 
its direction, an interesting future issue is the incorporation 
of wind information into our segmentation algorithms. By 
adopting anisotropic MRF in the prior, we intend to reflect this 
directional dependency of the clique potentials. In the practice, 
the smoothness parameter /3 is no longer considered to be a 
constant value but is clique-dependent according to the wind 
direction and/or velocity. For estimating the wind profile from 
the SAR data, state-of-the-art algorithms are used. We have run 
promissing simulations, are currently testing this extension on 
real data, and expect to obtain interesting results on a near 
future. 



Appendix 
Fitting a Mixture of Gamma Densities 

As stated in the main text of this work, when segmenting 
real SAR images, the adopted data model for each class is a 
finite Gamma mixture. For completely defining the data model 
for each class (for lightness of the notation we have dropped 
the class index), we need to estimate 2K Gamma parameters, 
= (6>i, . . . , Ok), with Os = (ag, Ag), and K a priori proba- 
bilities a = (ai, . . . , ax)). We infer (/) = (a, 0) by computing 
its ML estimate from a training set. The ML estimate is 
computed via an EM algorithm [?]. Fitting a Gamma mixture 
is addressed in [?]. However, the authors consider A/^-look 
SAR images meaning that the underlying random variables 
are the average of N independent and identically distributed 
exponential random variables, thus, having a Gamma density 
but with just one parameter free; if the mean is /i, then the 
variance is given by ji^ /N. We estimate both the mean and the 
variance for each Gamma distribution in the mixture, rendering 
the algorithm more adaptable to real measurements. 

The key point in the EM technique is the introduction of the 
so called missing data z, such that p(y|(/>) = J p {y ^ z\(j)) dz 
and p{y^z\(j)) is easier to manipulate than p{y\(j)). In the 
particular case of a mixture of densities, we will use as 
missing data a random variable, Zi, per site, with distribution 
p{zi = s) = as. It is interpretable as the probability of 
the s — th Gamma mode is selected at pixel i. The EM 
algorithm alternates between two steps: the E-step computes 
the conditional expectation of the logarithm of the complete 
a posteriori probability function, with respect to the missing 
variables, based on the actual parameter value. The M-step 
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updates the values of the parameters, by maximizing the ex- 
pression obtained in the E-step with respect to each parameter 
on turn, i.e., 

E-step : Q{^;^')=E{\ogp{y,z\^)\y,^'} (29) 
M-step : (/)^+^ = arg max Q (</>; (/)^) (30) 
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Denoting 



P{zi = s\yi,(l)*) 



(31) 



and taking into account that XI <^s - 1' then ^^+^ can be found 
among the stationary points of the Lagrangean 

N K / K \ 

with 

Ls^{Os,as) = log(A^O-log[r(a,)] 

+ log(7/,"^-^) - A,7/, + log(a,), (33) 

where A denotes a Lagrange multiplier. The expression for wl^ 
(see [?]) is given by 



w. 



(^\p{yM) 



^K 



Er=l<^(^i| ^r) 



(34) 



In the M-step, after differentiating C in order to the unknown 
parameters and setting the derivatives to zero, we obtain a 
closed solution for the updating of the a priori probabilities 
Qf^'s, but numerical iteration is needed for determining param- 
eters a^'s and A^'s of the Gamma densities. Expression ([35]) 
gives the update expression for a^'s. 



v^+l - 



1 



N 






(35) 



Equation (l36l) and (|37]) give the update expressions for the 
parameters A^'s and a^'s. 



A*+i = 






where 



\t\srN , 



N 



log iK) Ei=i wji + Ei=i log JyiWs 
Ei=i Ki 

Tias)' 



*K) 



(36) 

1 

(37) 
(38) 



is the psi function. We refer to the Appendix B of [?] for a 
very fast Newton procedure to compute the inverse of the psi(-) 
function. Expressions (l36l) and (|37]) are iteratively recomputed 
until convergence is obtained, starting from initial values 
computed from the observed data y. The initial parameter 
values are calculated in such a way, that the initial probability 
function is a sum of equidistant Gammas that span the most 
representative data range. The EM scheme converges in a few 
tens of iterations. 
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